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It is expected that the interplay between non-trivial band topology and strong electron correlation 
will lead to very rich physics. Thus a controlled study of the competition between topology and cor¬ 
relation is of great interest. Here, employing large-scale quantum Monte Carlo (QMC) simulations, 
we provide a concrete example of the Kane-Mele-Hubbard (KMH) model on an AA stacking bilayer 
honeycomb lattice with inter-layer antiferromagnetic interaction. Our simulation identified several 
different phases: a quantum spin-Hall insulator (QSH), a xy-plane antiferromagnetic Mott insulator 
(xy-AFM) and an inter-layer dimer-singlet insulator (dimer-singlet). Most importantly, a bona fide 
topological phase transition between the QSH and the dimer-singlet insulators, purely driven by the 
inter-layer antiferromagnetic interaction is found. At the transition, the spin and charge gap of the 
system close while the single-particle excitations remain gapped, which means that this transition 
has no mean field analogue and it can be viewed as a transition between bosonic SPT states. At one 
special point, this transition is described by a (2 + l)d 0(4) nonlinear sigma model (NLSM) with 
exaet SO(4) symmetry, and a topological term at exaetly © = tt. Relevance of this work towards 
more general interacting SPT states is discussed. 

PACS numbers: 71.10.-w, 71.10.Fd, 71.27.-ha 

I. INTRODUCTION 

The interplay between non-trivial band topology 
and strong electron interaction is expected to lead 
to a plethora of new physical phenomena in strongly 
correlated systems. Many exotic phenomena of in¬ 
teracting topological insulators (TI) have been pre¬ 
dicted/discovered, such as topological Kondo insula¬ 
tor^”^, fractionalized TI^’^, interaction-reduced classifi¬ 
cation of and interaction-driven anomalous topo¬ 

logical order at the boundary of TIs^^”^^. Besides 
fermionic systems, it was also proposed that bosonic 
systems can also form exotic states that are similar to 
fermionic TIs^^’^^, which are generally called the sym¬ 
metry protected topological (SPT) states. Unlike their 
fermionic counterparts, bosonic SPT states can only ex¬ 
ist in strongly interacting boson systems, and the inter¬ 
action must be carefully designed to avoid the ordinary 
superfluid and Mott insulator phases. These studies have 
tremendously broadened our understanding of quantum 
disordered states of matter and revealed the fundamental 
role topology plays in condensed matter systems. 

Quantum phase transition between different stable 
quantum disordered phases is another important subject, 
and in general it can be very different from the stan¬ 
dard Ginzburg-Landau (GL) phase transition paradigm. 

For example, one expects a phase transition between a 
(2 + l)d topological ordered state (Z 2 spin liquid^^) and 
an conventionally ordered phase (superfluid) is beyond 
the GL paradigm, and the Landau order parameter will 
acquire an enormous anomalous dimension. This phe¬ 
nomenon is confirmed by unbiased quantum Monte Garlo 
simulations^^’^^. In the non-interacting limit, the quan¬ 


tum critical point between two different topological insu¬ 
lators is usually described by a gapless Dirac/Majorana 
fermion, but the role of strong interaction at this tran¬ 
sition has not been fully explored, although we under¬ 
stand that in some particular cases interaction can gap 
out this quantum critical point and lead to a continuous 
curve connecting the two sides of the phase diagram^’^. 
Quantum phase transitions between bosonic SPT states 
were even less studied, and it was pointed out that most 
generally two bosonic SPT states can be separated by an 
intermediate phase^^’^^. 

With this in mind, it will be of great interest to inves¬ 
tigate a concrete example where in a strongly correlated 
fermionic SPT setup there is a purely interaction-driven 
phase transition between a topological insulator and a 
quantum disordered phase. Such a bona fide interaction- 
driven topological phase transition will have no mean- 
field (non-interacting) correspondence and provide the 
precious example of a controlled study of the interplay 
between non-trivial band topology and strong electron 
interaction. And this is what we will focus on in this 
paper. 

Here, we provide a concrete simple interacting fermion 
model that is studied by large-scale unbiased QMG sim¬ 
ulations. The results of this investigation provide us 
with the following desired phenomena: A bona fide 
interaction-driven quantum phase transition between 
topological insulator and a strongly interacting Mott in¬ 
sulator (a quantum disordered phase). We find that 
this quantum critical point is fundamentally different 
from the Tl-to-trivial quantum phase transition in the 
free fermion limit, in the sense that the fermions never 
close their gap at the transition, but emergent collective 
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bosonic degrees of freedom become critical. Thus we can 
view this transition as a transition between a bosonic 
SPT state and a trivial bosonic Mott insulator. And we 
demonstrate that at one special point, this transition is 
described by a (2 + l)d 0(4) nonlinear sigma model with 
exact SO {4) symmetry, and a topological term at exactly 
0 = TT. Moreover, we also employ the strange correlator 
proposed by Ref.^^ and tested in Ref.^^“^^ to diagnose 
the topological nature of the interaction-driven quantum 
phase transition between topological insulator and the 
strongly interacting Mott insulator. 


II. MODEL AND NUMERICAL METHOD 

A. AA-stacked bilayer KMH model with 
inter-layer AFM coupling 

In this work, we employ large-scale QMC simulations 
to investigate the AA-stacked bilayer KMH model with 
inter-layer AFM coupling, the Hamiltonian is given by, 
H = Htb + Hu + Hj^ as 

H =—t 

- 1 )^ 

+ ^ E - ^1,2^ - + Dl^2if] , (1) 

with Dii^ 2 i = Zlcr P denote the spin species 

t and I and ^ = 1,2 stand for the layer index in the 
AA-stacked bilayer system, as shown in Fig. 1. Htb de¬ 
scribes the tight-binding part of the Hamiltonian, includ¬ 
ing the nearest-neighbor hopping and the spin-orbit cou- 
phng36,37 Q^] 2 d the factor Vij = —Vji = ±1 depends 

on the orientation of the two nearest neighbor bonds that 
the electron traverses in going from site j to i, as shown 
in Fig. 1 (a). The in the spin-orbit coupling term 
furthermore distinguishes the t and | spin states with 
the opposite next-nearest-neighbor hopping amplitude. 
Throughout this work, we take t as unit of energy. The 
second term Hu describes the on-site Coulomb repulsion 
between electrons, and n^i = The electron fill¬ 

ing is fixed at half-filled, i.e., one electron per site on 
average. The third term Hj stands for the inter-layer 
antiferromagnetic spin interaction. As explained in de¬ 
tails in Appendix A, it is a faithful approximation of the 
full Heisenberg interaction ’ ^ 2 i- 

The Kane-Mele (KM) model preserves time-reversal 
symmetry Zj and its ground state is a quantum spin- 
Hall insulator with counter propagating edge states^^’^^. 
On the AA-stacked bilayer honeycomb lattice, the ground 
state of KM model is still a QSH insulator but with 



FIG. 1. (color online) (a) Illustration of AA-stacked honey¬ 
comb lattice and bilayer KMH model with inter-layer anti¬ 
ferromagnetic exchange interaction. The four-site unit cell is 
presented as the shaded rectangle. The gray and black lines 
indicates the nearest-neighbor hopping t on layer 1 and 2, 
respectively. The spin-orbital coupling term A, for one spin 
flavor, is shown by the red lines and arrows with I'ij — +1. 
The on-site Coulomb repulsion and inter-layer AFM coupling 
are represented by the shaded circle and rectangle, respec¬ 
tively. (b) Illustration of the xy-AFM Mott insulator phase, 
(c) Illustration of the inter-layer dimer-singlet phase. Shaded 
ellipses are the inter-layer spin singlets. 


two sets of counter-propagating edge modes. As for the 
symmetry, the model Hamiltonian in Eq. (1) has charge 
1/(1) X 1/(1) symmetry, which corresponds to charge con¬ 
servation on each individual layer. The SU{2) spin- 
rotational symmetry is broken down to U{1) by the spin- 
orbit coupling term, the residual 1/(1) spin symmetry cor¬ 
responds to the spin rotation in the xy plane. Therefore, 
most generally the total symmetry of the AA-stacked 
bilayer model is //(l)spin ^ [^(^) ^ ^(I)]charge ^ ^2 ’ 
which results in a Z classification. This is because in 
the noninteracting limit we can define a Chern number 
for spin-up and spin-down electrons separately, and time- 
reversal symmetry guarantees that these two Chern num¬ 
bers must be equal. Thus eventually the whole system is 
characterized by one Chern number, which can take arbi¬ 
trary integer values. When including interaction terms, 
we found that at the limit oiU = 0, Eq. (1) has a much 
higher SO (4) symmetry, which we will analyze in detail 
in Sec. HIB. 

With interactions, the KMH model on the monolayer 
honeycomb lattice has been studied by the Hartree- 
Eock mean-field theory^^, cluster (dynamic) mean-field 
theory^^”^^ as well as determinantal QMC simula¬ 
tions^^ Eor the bilayer model in Eq. (1), at the 
U = J = 0 limit, the system is a QSH insulator with spin 
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Chern number Cg = {C^ — C^)/2 = 2 where = +2 
and C]^ = —2 are the Chern numbers for spin-up and 
spin-down parts. In presence of finite interactions, i.e., 
the U — J phase diagram, one can expect that in the 
large U limit, the bilayer system will be driven from 
the QSH state into a xy-AFM ordered Mott insulator 
phase, through a continuous phase transition, similar to 
the KMH model on monolayer honeycomb lattice^^’^^“^^, 
and the phase transition should belong to the (2 -h l)d 
XY universality class^^’^^. At large J limit, the bilayer 
system should enter the inter-layer dimer-singlet phase 
with spin singlets formed on the inter-layer bonds due to 
strong antiferromagnetic coupling J. In the J ^ oo limit, 
the inter-layer dimer-singlet phase is a product state of 
the inter-layer singlets in which all the symmetries are 
preserved^^. Combining the spin-orbit coupling term A, 
on-site Coulomb repulsion U and inter-layer coupling J, 
one can expect very interesting competition occurring 
among the QSH, xy-AFM and inter-layer dimer-singlet 
phases, and it is the quantum phase transitions between 
these phases (some of which is of exotic topological na¬ 
ture) that we engaged great effort to unravel in this paper 
with unbiased large-scale QMC simulations. 


B. Projector quantum Monte Carlo method 

Projector QMC (PQMC) method is the zero- 
temperature version of determinantal QMC algorithm^^. 
PQMC method obtains the ground-state expectation val¬ 
ues of physical quantities by carrying out an imaginary- 
time evolution of some trial wavefunction, which is not 
orthogonal to the true many-body ground state. The 
ground-state expectation value of physical observable is 
calculated as follows. 


To determine the phase diagram for the bilayer model 
in Eq. ( 1 ), we first measure static physical quantities, 
such as the expectation values of energy densities (both 
total and each individual term in the Hamiltonian), dou¬ 
ble occupancy, and spin-spin correlation function. The 
x^-plane AFM order is expected to have ordering vec¬ 
tor r = ( 0 , 0 )^^’^^“^^, the transverse magnetic structure 
factor at P point is measured as, 

= iV (3) 

b'7 

where = 1 , 2 ,-- - ,Y run over all unit cells and 7 = 
1, 2, 3,4 stands for the four sublattices inside a unit cell. 
The staggered magnetic moment ms can be evaluated as 
ms = ^/S^y{T)/N. 

Next, to have the dynamical information of the system, 
such as the excitation gaps in single- and two-particle 
channels, we need to measure the imaginary-time single¬ 
particle Green’s function, 

= iV Y (4) 

ij^a 

where 7 is again the sublattice index and a is the elec¬ 
tron spin, and the imaginary-time spin-spin correlation 
function at F point, 

S^^r, r) = T (5) 

ijj 

and the imaginary-time inter-layer pair-pair correlation 
function in the charge channel 

+ At,(r)A,,) (6) 

ijS 


{ 6 )= lim 

0^+00 

where I^Pt) is the trial wave function and 0 is projection 
parameter. In all the simulations, to ensure that the al¬ 
gorithm arrives at the truly converged ground state of 
finite size systems, we choose 0 = 60/t,Ar = 0.05/t, 
in which Ar is the finite imaginary-time step applied in 
the Trotter decomposition of partition function. Dur¬ 
ing the simulations, we adopt the Hubbard-Stratonovich 
(HS) transformation with four-component Ising fields to 
decouple the interaction terms^^. Due to the fact that the 
two terms in Hj interaction do not commute, the system¬ 
atic error for all physical observables is at the order of 
0(Ar) (Trotter Error). During the simulation, we make 
sure the value of Ar is small enough and the QMC sam¬ 
pling of physical observables is large enough such that the 
results are numerically exact within well-characterized 
statistical errors. We have simulated different linear sys¬ 
tem size L = 3,6,9,12 (L = 15 for strange correlator), 
with A = as number of unit cells to extrapolate phys¬ 
ical observables to the thermodynamic limit. 


where Ai^ = is the inter¬ 

layer Cooper pair operator, it is defined on the two inter- 
layer bonds 6 = 1,2 of each unit cell i. At the r ^ 
oo limit, we access the asymptotic behavior G(k, r) oc 
g-Asp(k)r^ S^y(F^r) oc and P(r,r) oc in 

which Asp(k) is the single-particle excitation gap and 
A^, /Ac are the two-particle excitation gaps in the spin 
and charge channels for the interacting system^^. In our 
bilayer system, the minimum value of single-particle gap 
appears either at k = K or k = M depending on the 
parameters U and J, and we measure the spin and charge 
gaps at r point as it is the ordered wave vector for the 
gapless Goldstone modes. 

To diagnose the topological nature of the quantum 
phase transition, we employ the recently developed 
strange correlator method^^”^^. In the single-particle and 
two-particle (spin) channel, the correlation functions con¬ 
structed as. 


( 7 ) 
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where 4^^ = = 

A (iiiteger i as unit cell index), with 
k inside the BZ region, A, B are the sublattices in a 
unit cell in one layer and ^ the layer index, as shown 
in Fig. 1. The basic idea of the strange correlator is 
that, on the left hand side of the correlation function, the 
wave function \Q) is a trivial band insulator (with spin 
Chern number Cs = 0); on the right hand side of the 
correlation function, the projection operator guar¬ 

antees |T) = e“®^|TT) is the many-body ground state 
wave function of bilayer KMH Hamiltonian at certain J 
and U. If |T) is topologically nontrivial QSH state, i.e., 
there exit gapless edge modes at the spatial boundary of 
|T), then after a space-time rotation, will develop 

a singularity at certain symmetric momentum point k^: 
Ck ^ l/|k — ks|^, with a = 1 for noninteracting sys¬ 
tem, a < 1 otherwise. Based on the effective Lorentz 
invariant description of topological insulators^^, the 2D 
strange-correlator should behave very similarly to 

the (1 + l)d correlation functions at the boundary, en¬ 
dowed with a Luttinger liquid description in the presence 
of interaction. If |T) is on the other hand a topological 
trivial insulator, then the divergence in is no longer 

present as there is no single-particle edge modes on the 
boundary of |T). What’s more, the spin strange cor¬ 
relator ^SkAA different behaviors, depending on 

whether the gapless two-particle edge modes is present 
or not. For QSH insulator, should possess a di¬ 

verging behavior faster than ^ InL (the case in nonin¬ 
teracting system) with increasing system size L, while it 
should saturate to finite value (slower than ^ InL be¬ 
havior) in a topological trivial insulator. Thus, one can 
readily detect the topological phase transition in the sys¬ 
tem by monitoring the behavior of and 5'^^^. The 

strange correlation has been successfully applied in the 
QMC investigation of the topological phase transitions 
in the monolayer KMH model, the readers are referred 
to Ref.^^ for more details in its physical meaning and 
technical implementation. 


III. NUMERICAL RESULTS AND 
DISCUSSIONS 

A. Phase diagram 

The U — J phase diagram for A = 0.2t, 0.3t is shown in 
Fig. 2, and this is one of the main results of the paper. 
QSH, x^-AFM and inter-layer dimer-singlet phases are 
found from QMC simulations. Since there is only a net 
shift in the phase boundaries between A = 0.2t and 0.3t 
cases, we will focus on the detailed results for the A = 0.2t 
case in the following. The orange dotted line in Fig. 2 
denotes the J = 2U path which is studied in Ref.^^, we 
note that with more careful finite size scaling in this work, 
we found it actually goes through an intermediate AFM 
region. 


Three featuring observations about this phase diagram 
are in order. 



FIG. 2. (color online) U-J phase diagram for the AA-stacked 
bilayer KMH model with inter-layer antiferromagnetic cou¬ 
pling. Showing here are the phase diagram for A = 0.2t and 
A = 0.3t cases. Solid lines (violet, green and black) are the 
phase boundaries for the A = 0.2t case. The red solid dot at 
(Jc, U = 0) and red open dots at U = 0.25 and 0.5 and the 
green line goes through them highlights the interaction-driven 
topological phase transition between QSH and the dimer- 
singlet insulator phase. The orange dotted line highlights 
the J = 2U path which is studied in Ref. 50, it actually goes 
through a small AFM region. 

First of all, at small U {U < 0.5t for A = 0.2t) 
there is a direct phase transition from the QSH insu¬ 
lator to inter-layer dimer-singlet insulator (see details in 
Appendix D). Notice that since neither the QSH nor the 
dimer-singlet phase has symmetry breaking, all the sym¬ 
metries (such spin-rotation, charge conservation, time- 
reversal, and lattice symmetry, etc.) in the model Hamil¬ 
tonian Eq. (1) are preserved across this phase transi¬ 
tion, rendering it a bona fide topological phase transition 
driven purely by the inter-layer antiferromagnetic inter¬ 
action J. This is a very unique case and very different 
from the transitions in (interacting) topological insula¬ 
tors that have been studied before^^’^^“^^’^^, where the 
transitions are either driven by hopping parameters at 
free-fermion level^^’^^’^^’^^’^^, or after the transition the 
symmetry that protects the non-trivial band topology 
has been destroyed by interactions^^’^^“^^’^^. The nature 
of this exotic transition will be further discussed in the 
next section. 

Secondly, the region of xy-AFM phase is greatly ex¬ 
tended by an interesting collaboration between the on¬ 
site Coulomb repulsion U and the inter-layer AFM cou¬ 
pling J. At J = 0, for A = 0.2t, the QSH to xy-AFM 
phase transition occurs at U ~ 5.6(2)t^^, but as J in¬ 
creases, the phase boundary between QSH and xy-AFM 
moves towards smaller U, which means J and U both 
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prefer the AFM state, until J dominates over [/, after 
which the dimer-singlet phase takes over. The same phe¬ 
nomena is also observed for A = 0.31 case. 

Thirdly, for the direct phase transition from QSH 
phase to inter-layer dimer-singlet phase, we have ob¬ 
served signatures of continuous phase transitions. This 
can be seen from the inter-layer spin-spin correlation 
function per bond, shown in Fig. 3 (a) for L = 6 system 
with A = 0.2t (data with larger system sizes are shown 
in Appendix D). For various U values, as a function of 
J, the spin-spin correlation function changes from 0 to 
—3/4, with the latter signifying the formation of spin- 
singlet on every inter-layer bond. Moreover, according 
to the Hellmann-Feynman theorem, the spin-spin cor¬ 
relation function per bond is the first-order derivative of 
the total energy density over J. Combining the results of 
{Sii • S 2 i) presented in Fig. 3 (a) and in Appendix D, the 
continuous changing of the first-order derivative of the 
total energy density, with increasing J, suggests that the 
topological phase transition from QSH to dimer-singlet 
insulator phase is continuous (at least for U = 0). 

To further elaborate upon this point. Fig. 3 (b) and (c) 
show the first-order derivatives of expectation values of 
{Hj) per bond and {Hu) per site, over the parameter J. 
The peaks in Fig. 3 (b) indicate the QSH to dimer-singlet 
{U = 0) and xy-AFM to dimer-singlet (when U > 2t) 
phase transitions. The peaks in Fig. 3 (c) indicate not 
only the same transitions in Fig. 3 (b) at large J and 
small t/, but more interestingly, also the QSH to xy-AFM 
phase transitions at small J and large U (for U > 3t), 
as there are two peaks in the curves for U = 3t, 4t, 5t. 
The finite-size effects in the energy density derivatives 
are small, we only observe a slight shift of the phase 
boundaries for L = 9,12 systems, comparing with those 
for L = 6 system shown here. In the next section, we 
will present the finite-size scaling of the QMC results 
of the magnetic order parameter as well as the single¬ 
particle and spin excitation gaps across the topological 
phase transition between QSH and dimer-singlet. As we 
will see, the results hence obtained are consistent with 
those in Fig. 2 and Fig. 3 in this section. 


B. Topological phase transition 

1. Excitation gaps 

As mentioned in the preceding section (Sec. HI A), one 
of the most exciting features in the phase diagram (Fig. 2) 
is the exotic topological phase transition purely driven by 
the inter-layer antiferromagnetic interaction J, between 
the QSH and dimer-singlet phases. 

In a free-fermion system, topological phase transitions 
between SPT states are driven by tight-binding param¬ 
eters. The single-particle excitation gap will close to 
zero and reopen continuously at the transition, as long 
as the symmetries protecting the topologically nontriv¬ 
ial phase are still preserved. However, the topological 



j/t 


FIG. 3. (color online) (a) The inter-layer spin-spin correlation 
function for L = 6, A = 0.2t system with various U values, as 
a function of J. The continuous variation of this correlation 
function indicates the topological phase transition from QSH 
to dimer-singlet is a continuous one. At large J, the cor¬ 
relation saturates at —3/4 which signifies the formation of 
inter-layer dimer singlets, (b) First-order derivative of (Hj) 
per bond over J for L — 6, A = 0.2t. The peak in every curve 
explicitly indicates phase transition from QSH insulator (or 
xy-AFM) phase to inter-layer dimer-singlet phase, (c) First- 
order derivative of {Hu) per site over J. The peaks in these 
curves indicate all three possible phase transitions: QSH to 
dimer-singlet, QSH to xy-AFM and xy-AFM to dimer-singlet 
transitions. For U < 2t, the peak in d{Hu)ldJ corresponds 
to the QSH to dimer-singlet transition; for U > 3t, the two in¬ 
dependent peaks as a function of J correspond to the QSH to 
xy-AFM transition at small J and xy-AFM to dimer-singlet 
transition at large J. 
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FIG. 4. (color online) (a) Single-particle gap Asp(K) of A = 
0.2t,U = 0 as a function J. The inset shows the Asp(K) 
in J G [3.4t, 3.8t] region. We have checked that K point is 
indeed the minimum of single-particle gap in the whole BZ. 
As a function of J, the single-particle gap only shows a gentle 
dip near the topological phase transition, (b) Spin gap As of 
A = 0.2t, = 0 with increasing J. The inset is the spin gaps 

in J G [3.4t, 3.8t] region. The spin gap drops very fast and 
closes at the topological phase transition point Jc — 3.73(l)t. 

phase transitions in interacting systems seems to be much 
more complicated. Of course, they can still be driven 
by some tight-binding hopping parameter in the model 
Hamiltonian, such as the third-nearest-neighbor hop- 
ping4i,46,48,49, dimerized nearest-neighbor hopping^^“^^, 
Rashba spin-orbit coupling^^, and Kekule distortion^^ 
in the monolayer KMH model. In these cases, single¬ 
particle gap closes and reopens at the topological phase 
transition, just as their non-interacting counterparts. 
But, they can furthermore be driven purely by inter¬ 
actions, such as on-site Coulomb repulsion in mono- 
layer KMH inter-layer AFM ex¬ 

change coupling in AA-stacked bilayer KMH model^^, 
and more complicated form of interaction in interacting 
BHZ model^^ 

For interaction-driven topological phase transitions, 
the on-site Coulomb repulsion in the monolayer KMH 
model drives the QSH phase into an antiferromagneti- 
cally ordered phase with broken time-reversal and spin 



FIG. 5. (color online) Spin gap in J G [3.7t, 3.8t] region for 
A = 0.2t, H = 0 with L = 3,6,9,12 and the extrapolation 
by third-order polynomial. The inset shows the extrapolated 
spin gap as a function of J. 


rotational symmetries. Precisely speaking, this is still 
not the topological phase transition we are after in this 
paper: what we found here is an purely interaction- 
driven topological phase transition without any symme¬ 
try breaking on either side of the transition. Examples 
of this type of phase transition has been discussed in 1- 
dimensionaP^ and 2-dimensional^^ interacting systems. 
In Ref.^^, the single-particle gap remains gapped at the 
transition and it is the spin excitation gap and Cooper 
pair gap that close and reopen. This implies that in 
the low-energy limit such topological phase transition 
only involves bosonic degrees of freedom, allowing the 
fermionic excitations to be integrated out from the field 
theory 

Following Ref.^^, we perform a detailed study on the 
topological phase transition between QSH and dimer- 
singlet phases in the phase diagram of Fig. 2. To char¬ 
acterize this phase transition, we measured the single¬ 
particle gap, two-particle spin and charge gaps, as well 
as the strange correlator^^^^^ in the QMC simulations. 

The results of single-particle and spin gaps with in¬ 
creasing J are shown in Fig. 4 for A = 0.2t, 7/ = 0 in 
L = 3,6,9,12 bilayer systems. The raw data of the 
single-particle Green’s function and dynamic spin-spin 
correlation function are shown in Appendix B, the data 
are of very good quality, upon which we extracted the 
excitation gaps reliably. At 7/ = 0, the topological phase 
transition point is at Jc — 3.71 — 3.St in the phase dia¬ 
gram in Fig. 2. As shown in Fig. 4 (a), the single-particle 
gap only exhibits a very gentle dip around the topological 
phase transition point, which suggests the single-particle 
gap of the system remains open as a function of J. In 
contrast, we observe in Fig. 4 (b) that the spin gap de¬ 
creases rapidly in the vicinity of Jc as a function of system 
size L. The inset of Fig. 4 (b) shows the gap values in 
the region J 3.4t —3.8t. Within an even smaller region 
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of J G [3.7t,3.8t], we extrapolate the spin gap values for 
1/ = 3, 6 , 9,12 systems in l/L to estimate the spin gap in 
the thermodynamic limit, which is shown in Fig. 5. The 
main panel and inset of Fig. 5 deliver a clear message 
that the spin gap closing point is around Jc = 3.73(l)t. 
Furthermore, at = 0 as a function of J, we don’t find 
a stepping of xy-AFM order by finite size extrapolation 
of the transverse magnetic structure factor (see details 
in Appendix D). At J = 3.8t, the spin gap values for 
L = 9 system and L = 12 are almost the same, indi¬ 
cating that the thermodynamic limit is already reached 
and the spin excitations are well gapped here (spin-spin 
correlation in real space is exponentially short-ranged). 
After the topological phase transition, the bilayer sys¬ 
tem enters the inter-layer dimer insulator phase, which 
is schematically shown in Fig. 1 (c). 

Combining the results for single-particle and spin gaps, 
we find that the topological phase transition driven by 
inter-layer AFM coupling in our bilayer system is funda¬ 
mentally different from those controlled by the hopping 
parameters, with and without interactions^^’^^^^^’^^. 

Also, at U = 0 , we have observed that the charge gap 
Ac and spin gap As are numerically identical (with dif¬ 
ference only up to O.OOlt). There is actually a deep theo¬ 
retical reason for the equality between these two-particle 
gaps: it is due to an exact 5'0(4) symmetry at [/ = 0 (see 
Appendix C), which rotates the xy-AFM (spin) fluctu¬ 
ation Ni = ^(—and the pairing (charge) 
fluctuation = ciiia'^C 2 i like an 0(4) vector: 

m = (7Vf,ImA,,ReA,,7Vf). ( 8 ) 

Therefore both the spin and the charge excitation gaps 
close identically at the transition point. To better under¬ 
stand the SO (4) symmetry, we may define two fermion 
doublets fia (a =t,i): 


M = 


ciit A 


J 


(9) 


where we have Xcr = (~ 1)^5 and tij = t for hoppings 
on the NN bonds and tij = iX for SOC on the next- 
nearest-neighbor (NNN) bonds. Under arbitrary SU{2) 
rotation of operator as i-G the Pi opera¬ 

tor in Eq. (12) is invariant since we have fiP^‘^{fiaP 

, and the equality for 

2 x 2 SU{2) matrix 1/^- Besides, the hopping term 
flfjtijfjcr is explicitly invariant under the SU{2)a- ro¬ 
tation of U(j- Combining them, the Hamiltonian in 
Eq. ( 12 ) has independent SU{2)-^ and SU{2)^ symme¬ 
tries for spin up and down channels, respectively. Thus, 
the >50(4) SU{2)^ x SU{2)_i^ symmetry for the bilayer 

model in Eq. (1) under U = 9 condition, which can be 
expressed in Eq. (11), is explicit. 

Physically, the SO {4) symmetry rotates the four com¬ 
ponents of rii defined in Eq. ( 8 ) to one another. As a 
result, the x^-AEM order should be exactly degenerate 
with the inter-layer spin-singlet 5-wave superconducting 
order under U = 0 condition due to the SO {4) symmetry, 
which also indicates the identical excitation gaps corre¬ 
sponding to these two orders, i.e., spin gap and charge 
gap. 


2. Theoretical understanding 

In our phase diagram, the fermionic single-particle 
gap never closes with finite A, while the two-particle, 
collective, bosonic modes (spin and charge gaps) both 
close at the QSH-to-dimer-singlet phase transition, this 
means that at low energy this model can be well- 
approximated by a bosonic model. Indeed, Ref.^^ demon¬ 
strated that many bosonic SPT states can be constructed 
from fermionic topological insulators/superconductors by 
confining the fermionic degrees of freedom. In our case, 
we propose that the bosonic sector of our phase diagram, 
at U = 0 , can be described by the following nonlinear 
sigma model (NLSM) field theory^^: 


Then the 0(4) vector can be written as 

-i'r^)/4 + (10) 

where are the Pauli matrices acting on the /- 

fermion doublets. The SO(4) group is naturally factor¬ 
ized to SU{2)^ X SU{2)^ as right and left isoclinic ro¬ 
tations, under which the fermion transforms as i-G- 

fi^Ua with Ua G SU{2)cr for both a =t 5 i- The model 

Hamiltonian in Eq. (1) at U = 0 can be written in terms 
of the SU{2)^ X SU{2)i singlets as 

H = xMlUjfja + h.c.) - J +PiPhm 

i,j,cr i 

with 

Pi = li-^rT.flir\fl)'^ ( 12 ) 

a 


S = [ d^^xdr + ^iabcdn°'dxn'’dyn‘^drn‘^, (13) 

J 9 “3 

where U 3 = 27r^ is the volume of a three dimensional 
sphere with unit radius. We will focus on the phase with 
large g, namely the vector n is disordered. Eq. (13) is 
exactly the same field theory introduced by Ref.^^’^^ to 
describe 2d bosonic SPT states, and the physical meaning 
of the four component vector field n was given in Eq. ( 8 ). 
As we show explicitly in Appendix C, the model Eq. ( 1 ) 
at U = 0 has exactly >50(4) symmetry, thus we do not 
need to turn on any anisotropic term to Eq. (13). When 
we move away from the point U = 0 , an anisotropy needs 
to be turned on to split the degeneracy between ( 77 - 1 , 77 - 4 ) 
and ( 712 , 77 - 3 ). 

The phase diagram and renormalization group flow of 
the (1 + l)d analogue of Eq. (13) were calculated explic¬ 
itly in Refs.^^^^^; and it was demonstrated that the en¬ 
tire phase 0 < 0 < TT is controlled by the trivial fixed 
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point 0 = 0 , while the entire phase tt < 0 < 27 r will 
flow to the fixed point 0 = 27r. The phase diagram of 
Eq. (13) was studied in Ref.^^, and again in the disor¬ 
dered phases (phases with large 0 = tt is the quantum 
phase transition between the two phases with 0 < 0 < tt 
and TT < 0 < 27r, the stable fixed point 0 = 27r describes 
a bosonic SPT state in (2 + 

The physical meaning of the fixed point 0 = 27 r be¬ 
comes explicit when we create a vortex of A, i.e. the 
vortex of ( 712 ,^ 3 ), then this vortex will acquire spin -1 
due to the 0— term at 0 = 27r, which is consistent with 
two copies of quantum spin Hall insulator with con¬ 
servation. Also, at the fixed point 0 = 27 r, the boundary 
of Eq. (13) is a (l + l)d 0(4) NLSM with a Wess-Zumino- 
Witten term at level- 1^^’^^, whose SO (4) symmetry fac¬ 
torizes into SU{2)lX SU{2)r {SU(2) symmetries for left 
and right moving modes respectively), where SU ( 2 ) l and 
SU(2)r precisely correspond to SU{2)^ and SU{2)^ in¬ 
troduced in the previous subsection. Thus the field the¬ 
ory Eq. (13) does match with the all the desired physics 
of our lattice model. In a later paper by some of us^^, we 
demonstrate that the boundary state of our lattice model 
will be driven into a purely bosonic conformal field the¬ 
ory, in the sense that all the fermionic modes are gapped 
by interaction, but bosonic modes are gapless. And the 
remaining gapless bosonic modes at the boundary are 
precisely described by the boundary states of Eq. (13). 

In Eq. (13) 0 = tt is the quantum phase transition 
between the SPT and trivial phases, and in our phase 
diagram 0 = tt corresponds to the direct QSH-to-dimer- 
singlet phase transition. Thus our lattice model actually 
provides a way to simulate the topological field theory 
Eq. (13) in QMC without sign-problem. 


3. Strange correlator 

Let’s now turn to understand the topological phase 
transition from QSH to dimer-singlet phases from the 
perspective of edge states. At U = 0 , in the QSH phase 
with J < Jc^ there exists two pair of gapless edge modes 
on the boundary of the bilayer KMH system, i.e., the spin 
Chern number Cs = 2. When J > Jc, the system is the 
dimer-singlet state, it is a topologically trivial product 
state hence the edge states are no longer present, i.e., 
spin Chern number Cg = 0. Therefore, the change of the 
topological nature from QSH to dimer-singlet can be seen 
from the presence/absence of the gapless edge states. 

In the QMC simulations, one can explicitly probe 
the spatial edge by applying open boundary condition 
(OBC), but in interacting systems, OBC usually has very 
strong finite-size dependence. Moreover, to be able to see 
the edge mode, one further needs to analytically continue 
the imaginary time correlation functions to have the spec¬ 
tra in real-frequency, but it is well-known that analytical 
continuation usually generates ambiguous results to the 
fine features of the spectra. Hence, to avoid such diffi¬ 
culties, recently there is a new diagnosis dubbed strange 



L 

FIG. 6. (color online) (a) The inverse amplitude of single¬ 
particle strange correlator l/|Cj^yi^| along the high-symmetry 
path for various J = 3.0t, 4.0t. (b) The spin strange correlator 
S^aa various J values as a function of linear system size 

L. 


correlator, that has been proposed/tested successfully in 
probing the edge states from static, bulk wave functions 
with periodic boundary condition^^”^^. 

As explained in the Sec. HB, whether the gapless 
edge modes is present in the bilayer system or not can 
be signified by the divergence of the single-particle and 
spin strange correlator, which are shown in Eig. 6 for 
A = 0.2t, 1/ = 0. Erom the single-particle strange corre¬ 
lator results in Eig. 6 (a), for J = 3t {J < Jc)^ \^\^ab\ 
of the bilayer KHM mode is diverging at M point, corre¬ 
spondingly, ^/\CI^ab\ vanishes in a power-law (the expo¬ 
nent a is almost 1) to zero. The data point of ^/\CI^ab\ 
exactly at k = M is a finite-size effect due to the imple¬ 
mentation of strange correlator in QMC and has been ex¬ 
plained thoroughly in Ref.^^. But when J = At (J > Jc), 
the divergence of is removed, hence the 1 /|C^^^| 

is no longer vanishing at M point, resembling the single¬ 
particle edge modes in QSH being gapped out due to the 
inter-layer antiferromagnetic interaction J. As for the 
spin strange correlator shown in Eig. 6 (b), S^aa diverges 
with increasing L oA J < Jc^ which is faster than the InL 
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behavior at the noninteracting limit J = 0. These re¬ 
sults indict the existence of gapless, spin (bosonic) edge 
modes^^. On the contrary, simply saturate to finite 

values when J > Jc, suggesting the absence of gapless 
edge modes. Combining the results of the strange cor¬ 
relator in both single-particle and two-particle channels, 
the QSH phase (J < Jc) in the bilayer model has gap¬ 
less edge modes (bosonic), while they are absent in the 
dimer-singlet insulator phase, highlighting the topologi¬ 
cal phase transition. 

C. xy-AFNL order 

The xy-AFM order in the phase diagram Fig. 2 corre¬ 
sponds to the ordered phase g < gc in Eq. (13), with an 
extra anisotropy term that favors (ni,n 2 ) over ( 77 . 3 , 77 . 4 ). 



FIG. 7. (color online) (a) Finite-size extrapolation of the 
transverse magnetic structure factor for L = 3,6,9,12 sys¬ 
tems, the fits are third-order polynomial in 1/L. The param¬ 
eters are A = 0.2t, 1/ = 2t and J G [2.0t,2.7t]. Inset shows 
the extrapolated staggered magnetic moment ms as a func¬ 
tion of J. (b) The spin gap for L = 3, 6, 9,12 systems and its 
extrapolated thermodynamic limit (TDL) values for the same 
parameter set. 


In the phase diagram of Fig. 2, one finds the region 
of xy-AFM phase is greatly extended by an interesting 
collaboration between the on-site Coulomb repulsion U 
and the inter-layer AFM coupling J. Intuitively, the U 
term favors the xy-AFM state, while the J term favors 
the dimer-singlet state. With increasing 7/, the QSH to 
xy-AFM and xy-AFM to dimer-singlet phase transition 
points all move towards smaller J. This can be under¬ 
stood as following: the x^-AFM phase is triggered by the 
intra-layer antiferromagnetic coupling Jintra oc t^/7/, the 
dimer-singlet phase is triggered by the inter-layer J, their 
phase transition is determined by the ratio J/Jintra 5 since 
we get a smaller Jintra for larger 7/, the critical J for the 
phase transition to dimer-singlet is therefore reduced. 



1/L 



FIG. 8. (color online) (a) Finite-size extrapolation of the 
transverse magnetic structure factor for L — 3,6,9,12 sys¬ 
tems, the fits are third-order polynomial in 1/L. The param¬ 
eter sets are A = 0.2t, 7/ = 4t and J G [0.6t, 2.It]. Inset shows 
the extrapolated staggered magnetic moment ms- as a func¬ 
tion of J. (b) The spin gap for L = 3, 6, 9,12 systems and its 
extrapolated thermodynamic limit (TDL) values for the same 
parameter set. 

Let us be more quantitative about the phase bound¬ 
ary. For the monolayer KMH model with A = 0.2t, the 
system enters the xy-AFM phase at Uc = 5.6(2)t^^. In 
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the presence of inter-layer J, QMC results reveal that the 
xy-AFM phase can be well established even at ^ 2t. 
As shown in Fig. 7 (a), for the magnetic structure fac¬ 
tor for L = 3, 6, 9,12 systems and their extrapolation to 
thermodynamic limit in J G [2.0t,2.7t], the extrapolated 
S^y(T)/N takes nonzero values for J = 2.3t, 2.4t, 2.5t 
(see the inset of Fig. 7 (a)). To further confirm the long- 
range magnetic order, we have also measured the spin 
gap and the results are shown in Fig. 7 (b). The ex¬ 
trapolated spin gap at J = 2.3t, 2.4t, 2.5t are zero and 
corresponds to the Goldstone mode associated with the 
xy-AFM long-range order. Combining the data in Fig. 7 
(a) and (b), it’s very convincing that the long-range xy- 
plane magnetic order already appears at ^ 2t, almost 
3 times smaller than that of the J = 0 case. 

When the on-site Coulomb repulsion is further in¬ 
creased to U = 4t, at A = 0.2t and J G [0.6t, 2.It], there 
are two phase transitions (QSH to xy-AFM and xy-AFM 
to dimer-singlet) as J increases. These can be detected 
by measuring the magnetic structure factor and the spin 
gap as well, the results are show in Fig. 8. Fig. 8 (a) shows 
that the system is in xy-AFM phase in J G [l.Ot, 1.7t] by 
finite size extrapolation. The spin gap result in Fig. 8 (b) 
is well consistent with it, as the spin excitations are gap¬ 
less in the thermodynamic limit in J G [l.Ot, 1.7t]. When 
J < l.Ot, the system is inside the QSH insulator where 
the spin excitations are gapped, and when J > 1.7t, the 
system is inside the dimer-singlet phase where the spin 
excitations are gapped as well. 


IV. SUMMARY AND OUTLOOK 


of this transition. It seems the ordinary calculation tech¬ 
niques such as 1/N or e—expansion both fail here, be¬ 
cause Eq. (13) is defined solely for (2 -j- l)d and 0(4) vec¬ 
tor. How do we compare the critical scaling behavior of 
the spin gap measured in Fig. 5 to theoretical calcu¬ 
lations based on Eq. (13) is an interesting open question, 
which we will leave to future study. 
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In this work, we have found a bona fide interaction- 
driven quantum phase transition between topological in¬ 
sulator and a strongly interacting Mott insulator (dimer- 
singlet). This quantum critical point is fundamentally 
different from the Tl-to-trivial quantum phase transi¬ 
tion in the non-interacting limit, in the sense that the 
fermions never close their gap at the transition, instead, 
emergent collective bosonic degrees of freedom become 
critical. We also employ the strange correlator pro¬ 
posed/tested in Ref.^^“^^ to diagnose the topological na¬ 
ture of the quantum phase transitions. 

In principle the exotic topological phase transition we 
found in this paper can be generalized to all higher di¬ 
mensions. What we need to find is a higher dimension 
fermionic topological insulator/superconductor that can 
be mapped to a bosonic SPT state after confining the 
fermionic degrees of freedom, then in principle the simi¬ 
lar type of SPT-trivial phase transition with gapless bo¬ 
son modes but no gapless fermion mode can be found 
in these cases. A construction of these models in higher 
dimensions was discussed in Ref.^^. 

Although we have identified the field theory that de¬ 
scribes this interaction-driven direct Tl-to-trivial quan¬ 
tum phase transition in Eq. (13), we do not yet have a 
controlled analytical calculation for the universality class 


Appendix A: approximate Heisenberg interaction 


In Sec. HA, we mention the inter-layer antiferromag¬ 
netic interaction in our Hamiltonian is a faithful approx¬ 
imation of the full antiferromagnetic Heisenberg interac¬ 
tion. Here we elaborate more upon this point. 

The inter-layer interaction term Hj in Eq. (1) can be 
written as summation of the following term on all inter¬ 
layer bonds. 


Qi = \ 


{Du,2i - - {Du,2^ + lAifiif -(Al) 


There is an operator identity relates Qi with full Heisen¬ 
berg exchange coupling^^, it reads 


Sii • S2i = Qi — - (hi,i — l)(h2,i — 1) — 1 , (A2) 


so the difference between Qi and Sii • S 2 i is at the part 
[(f^i,i — l)(^ 2 ,i ~ 1) ~ 1] (where indexes 1,2 stand for 
layers and integer i for lattice site), but since our system 
is half-filled, the expectation value of (hi^i) = (h 2 ,i) = 1, 
i.e., the charge fluctuations are small. This term can be 
safely considered as a constant. 
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Moreover, it is easy to see Si^ • S 2 i and Qi share the 
same eigenstates and their eigenvalues are different only 
up to a 1/4 shift. The eigenstates for Si^ • S 2 i and Qi are 
spin singlet and three-fold degenerate spin triplet states, 

iV’0,+0) = (I t)ii 1)2 -1 i)ii 1)2) 

iV’1,+1) = (I t)ii 1)2) 

iV’1,+0) = (I t)ii 1)2 +1 i)ii 1)2) 

|^/’i,-i) = (U)i| 1)2) • (A3) 

for Sii • 822 , it’s well known that 

3 

Sli • S 22 |' 0 O,o) = -^|'0o,o) 

Sli • S2i\iJl,m) = m = 0, ±1. (A4) 

for Qi interaction, it’s simple to show 
Qi\i^O,o) = -1 • |'0o,o) 

Oil'll,m) = +0 • m = 0 , ± 1 . (A5) 


In terms of implementation in the PQMC simulations, 
for the Qi term, we can directly apply the following 
Hubbard-Stratonovich transformation to transform the 
Qi term into free fermion system coupled to 4-component 
Ising fields. 


exp 


J 




= I H + O [(Ar)^] , 


l=±l,±2 


exp 


J 


—(-Dli,2i + D\^ 2 iY 


= I II + O [(Ar)^] , (A6) 


Z=±l,±2 


with = a/AtJ/S. For the full 812-822 interaction term, 
we need to rewrite it into summation of Qi interaction 
term, the on-site attractive interaction (the second term 
in Eq. A7) and the inter-layer density-density attractive 
interaction (the third term Eq. A7) as follows, 

Sii ■ S2i = Qi - I [{ni,i - if + (n2,i - 1 )^] 

- +^2.i - 2)^. (A7) 

The problem here is that with J > 0 (the antiferro¬ 
magnetic interaction), the simultaneous presence of all 
these three terms will generate minus sign problem to the 
QMC simulation under 1/ > 0 condition as the model in 
Eq. ( 1 ), which effectively means that there is no way to 
perform QMC simulation with the full Heisenberg inter¬ 
action term for large systems. Though the QMC simula¬ 
tions applying the full J term as Eq. A7 for the bilayer 


model under U = 0 condition is free from sign problem, 
only keeping the Qi in 8 i^ • 822 interaction during the 
QMC simulations is still a good approximation, since the 
single-particle gap is always finite with U = 0 and arbi¬ 
trary J parameter. 


Appendix B: Raw data for dynamic correlation 
functions 

In Sec. HIB, we present the single-particle as well as 
the spin excitation gaps at the topological phase tran¬ 
sition between QSH and dimer-singlet phases. Here we 
show some raw data for imaginary-time single-particle 
Green’s function and spin-spin correlation function, to 
provide the evidence that the extrapolated excitation 
gaps are in good numerical quality. 



FIG. 9. (color online) Single-particle Green’s function for 
A = 0.2t, U = 0, J = 3.73t with L = 3, 6, 9,12 at K point in 
(a) linear scale and (b) in semi-logarithmic scale. 

Fig. 9 and Fig. 10 are the raw data of the single¬ 
particle Green’s function G(K, r) and the dynamic spin- 
spin correlation function 5'^^(r,r), with parameter set 
A = 0 . 2 t, U = 0, J = 3.73t. According to Fig. 4, this is 
exactly at J = Jc. In Fig. 9(a), we can observe the single¬ 
particle gap at K point decay very fast in imaginary time 
r. In Fig. 9 (b), with a semi-logarithmic scale, we can 
see the size of the single-particle gap almost converge 
to its thermodynamic limit value for L = 9,12 systems. 
Such fast decay and quick convergence with finite system 
size actually means the single-particle gap is indeed finite 
and large at the topological phase transition. In fact it 
is about 0.7t at the transition point. 

On the other hand, we can observe that the raw data 
for dynamic spin-spin correlation function in Fig. 10 (a) 
decay slower with r. And in Fig. 10 (b) with a semi- 
logarithmic scale, 5'^^(r,r) shows very good straight 
lines in imaginary time r and we can hence extract the 
spin gap value with very high accuracy. In fact, the 1/L 
finite size scaling of the spin gap at J = 3.73t gives rise 
a vanishing spin gap in the thermodynamic limit. 
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FIG. 10. (Color online) Dynamic spin-spin correlation func¬ 
tion for A = 0.2t, U = 0, J = 3.73t with L = 3, 6, 9,12 at T 
point in (a) linear scale and (b) in semi-logarithmic scale. 

Appendix C: The SO{4) symmetry 

As mentioned in Sec. II A, the bilayer KMH model 
given by Eq. (1) has the C(l)spinX [f 7 (l)x[/(l)]chargeXi^ 2 ’ 
symmetry in general. However when the model parame¬ 
ters are tuned to certain special combinations, the model 
can have larger symmetries. In this appendix, we will fo¬ 
cus on the various unitary symmetries of the model. The 
anti-unitary time-reversal symmetry Zj is always pre¬ 
sented and will be omitted in the following discussion. 

To understand the unitary symmetries systematically, 
let us first introduce three sets of competing orders (in 
terms of fermion bilinear operators): 

SDW: Ni = 

SC: Ai = cui<TJ'c 2 i, (Cl) 

Exciton: A = (-l)*ct'^ 2 i, 

where c^i = (c^^^, c^^)''’ is the fermion operator on site i 
of the ^ layer. ( — 1)^ and (—1)* respectively denote the 
staggered sign factors between the layers and between 
the sublattices. These competing orders anti-commute 
with each other, and can be organized into an 0(7) 
vector: Qi = (Af, Nf , Af, Re , Im , Re , Im ). 

Then one can introduce the SO{7) group on each site i 

that rotates the vector Qi. The generators of the SO{7) 
group are given by the following commutators (for a < b 
and a, 6 = 1, • • • ,7) 

= (C2) 

The fermion operator transforms under the SO(7) rota¬ 
tion (parameterized by Oab ^ as 

c^ia -)■ ex.p{i8abTf)c^ia expi-WabTf). (C3) 

The model Hamiltonian in Eq. (1) can not achieve this 
SO{7) symmetry, but its achievable unitary symmetries 


are all subgroups of this SO{7). Different choices of 
the model parameters breaks the SO{7) symmetry dif¬ 
ferently. 

To see how the SO(7) symmetry is broken explicitly 
by the Hamiltonian, we can calculate the commutator 
of the Hamiltonian H with the global SO{7) generators 

pa6 — pa6. 

= (C4) 

^ Q 

means the Hamiltonian has the symmetry that 
rotates and Q^. In general, is a linear combina¬ 
tion of operators with the model parameters X, U and 
J as coefficients: 

cab ^ ^cf + xcf + UCtI’ + JCf. (C5) 

cab ^ Cab^ C^b gj.g complicated operators whose 

detail expressions are not of much interest. We only need 
to extract the coefficients of linearly independent opera¬ 
tors, which are concluded in Tab. I. 

Most generally, only 3 (out of 21 ) SO{7) generators 
ri2, commute with the Hamiltonian, as = 

^45 _ (j67 _ They generate the 7/(l)spin x [^(1) x 
7/(1)]charge Symmetry group. However, when U = 0, we 
have = 0 in addition, which 

enlarges the symmetry group to 5'0(4) x 7/(1). The 
5'0(4) symmetry rotates the xy-SDW order and the SC 
pairing order as an 0(4) vector (A^, Re A, Im A), 

which involves particle-hole transformations. The 7/(1) 
symmetry rotates the exciton order (ReO,ImO) and 
corresponds to the conservation of the charge differ¬ 
ence between the layers. When J = 27/ 7 ^ 0, we 
have = 0 , which enlarges the symmetry 

group to SU(2)xU(l)spinXU(l)charge as mentioned in 
Ref.^^. When the interaction is completely turned off 
as 7/ = J = 0 , the model has 5'0(4) x SO(3) sym¬ 
metry. On the other hand, in the absence of the spin- 
orbital coupling, i.e. A = 0 , the model has even richer 
symmetry structures, as the spin SU{2) symmetry is re¬ 
stored. Under generic interaction, the symmetry group 
is 5'7/(2)spin X [7/(1) X 7 /(l)]charge 5 which can be enlarged 
to SO{b) X 7/(1) at 7/ = 0, or another SO{b) x 7/(1) at 
J = 27/, or 5'0(4) x SU{2) at J = 0. 

Appendix D: The topological phase transitions at 
small U region 

In Sec. HI, we present/discuss in detail the results 
about the J-driven topological phase transition without 
spontaneous symmetry breaking, including the energy 
derivatives, excitation gaps, strange correlator and quan¬ 
tum field theory correspondence. In this part, we show 
more numerical data about the topological phase transi¬ 
tion at small U region. 

As we have mentioned, the x^-AFM order is absent 
around the topological phase transition at small U re¬ 
gion. In Fig. 11, the extrapolation of structure factors of 
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J-2U 

SG 

Re A 




0 

J, A 

J, A 

ImA 





J, A 

J, A 

Exc. 

ReD 

ImD 



0 


TABLE L Linearly independent coefficients in For exam¬ 
ple, in row , column ReH, the number J — 2U and A mean 
that the commutation between H and ^ [N ^, ReH] contains 
an operator with coefficient J — 21/, and another operator 
with coefficient A. The details of the form of the operators 
are not shown. 



FIG. 11. (color online) Extrapolation of structure factor 
S^y{T)/N of xy-AFM order over 1/L for (a) 1/ = 0, (b) 
U = 0.251, (c) U = 0.501 and (d) U = 1.001 over inverse sys¬ 
tem size 1/L, at A = 0.21. The data points with error bars in 
the insets are the extrapolated values in the thermodynamic 
limit. 


xy-AFM order over 1/L for U = 0,0.251,0.501,1.01 are 
shown. From the results in Fig. 11, The xy-AFM order 
is explicitly absent for 1/ = 0 and U = 0.25t, correspond¬ 
ing of which the topological phase transition points are 
Jc/t = 3.731 and Jc/t = 3.541. For U = 0.51, only a 
single point of J/t = 3.37 has nonzero xy-AFM order 
applying the step size AJ = 0.011 during QMC simula¬ 
tions. Considering the numerical error existing in QMC 
simulations, it’s reasonable to terminate the xy-AFM or¬ 
dered phase at 1/ = 0.51 in the phase diagram presented 


in Fig. 2. For U = 1.01, the extrapolated xy-AFM order 
(inset of Fig. 11 (d)) is nonzero in 3.00 < J/t < 3.07 
region, which explicitly demonstrates the stepping in of 
xy-AFM ordered phase between the QSH insulator and 
inter-layer dimer-singlet insulator. Based on the results 
in Fig. 11 (a th- 
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FIG. 12. (color online) The inter-layer spin-spin correlation 
functions around the topological phase transitions for 1/ = 0 
with A = 0.21. 



J/t 


out spontaneous symmetry breaking is well established 
for 1/ = 0 and U = 0.251 case, i.e. finite U. 

Another question is whether the topological phase 
transition at small U is of first-order or continuous. Due 
to the fact that there is no nonzero local order parameter 
across the phase transition, to solve this problem thor¬ 
oughly is not easy. However, to resolve this problem as 
best as we can, we have measured the inter-layer spin- 
spin correlation function {Su • S 2 i) for 5 different system 
sizes at t/ = 0 as a function of J. As discussed in the 
main text, this quantity can be taken as the first-order 
derivative of ground state energy over J parameter of the 
model in Eq. (1). Depending on whether this quantity is 
continuous or not around the quantum phase transition 
in thermodynamic limit, we can determine the order of 
the transition. 

The results of {Su • S 2 i) for [/ = 0 across the topo¬ 
logical phase transition are shown in Fig. 12. We can 
observe that {Su • S 2 i) has almost reached the converged 
values already in L = 12 system, i.e. values at thermody¬ 
namic limit. This is rather reasonable since the fermionic 
channel of the system is always gapped and the finite-size 
effect should should not be so strong. Most importantly, 
we indeed observe that {Su'S 2 i) changes smoothly across 
the topological phase transition, which suggests a contin¬ 
uous phase transition. 
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